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Abstract 

We  consider  the  steady-state  equations  for  a  compressible  fluid.  For  low-speed  flow,  the 
system  is  stiff  because  the  ratio  of  the  convective  speed  to  the  speed  of  sound  is  quite  small. 
To  overcome  this  difficulty,  we  alter  the  time  evolution  of  the  equations  but  retain  the  same 
steady-state  analytic  equations.  To  achieve  high  numerical  resolution,  we  also  alter  the  artifi¬ 
cial  viscosity  of  the  numerical  scheme,  which  is  implemented  conveniently  by  using  other  sets 
of  variables  in  addition  to  the  conservative  variables.  We  investigate  the  effect  of  the  artifi¬ 
cial  dissipation  within  this  preconditioned  system.  We  consider  both  the  nonconservative  and 
conservative  formulations  for  artificial  viscosity  and  examine  their  effect  on  the  accuracy  and 
convergence  of  the  numerical  solutions.  The  numerical  results  for  viscous  three-dimensional  wing 
flows  and  two-dimensional  multi-element  airfoil  flows  indicate  that  efficient  multigrid  computa¬ 
tions  of  flows  with  arbitrarily  low  Mach  numbers  are  now  possible  with  only  minor  modifications 
to  existing  compressible  Navier-Stokes  codes.  The  conservative  formulation  for  artificial  viscos¬ 
ity,  coupled  with  the  preconditioning,  offers  a  viable  computational  fluid  dynamics  (CFD)  tool 
for  analyzing  problems  that  contain  both  incompressible  and  compressible  flow  regimes. 


*This  research  was  supported  in  part  by  the  National  Aeronautics  and  Space  Administration  under  NASA  Contract 
No.  NASl-19480  while  the  first  author  was  in  residence  at  the  Institute  for  Computer  Apphcations  in  Science  and 
Engineering  (ICASE),  NASA  Langley  Research  Center,  Hampton,  VA  23681-0001. 
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1  Introduction 

In  the  past  few  years,  several  preconditioning  methods  have  appeared  in  the  literature  [1-4]  with  the 
aim  of  solving  nearly  incompressible  flow  problems  with  numerical  algorithms  that  were  designed 
for  compressible  flows.  The  development  of  these  methods  are  motivated  by  two  main  observations. 
First,  flow  problems  exist  that  contain  both  compressible  and  incompressible  flows  simultaneously; 
that  is,  part  of  the  flow  region  can  be  considered  to  be  incompressible  with  locally  low  Mach 
numbers,  whereas  significant  compressibility  effects  occur  in  other  regions  of  the  flow.  A  typical 
example  in  aerodvnamics  is  the  flow  over  a  multielement  airfoil  near  maximum  lift.  Surface  heat 
transfer  or  volumetric  heat  addition  can  also  introduce  compressibility  effects  in  low-speed  flows. 
Second,  it  is  preferable  to  use  existing  compressible  flow  codes  over  the  broadest  range  of  flow 
conditions  possible  for  ease  of  use  and  consistency  reasons. 

The  difficulty  in  solving  the  compressible  equations  for  low  Mach  numbers  is  attributed  to  the 
large  disparity  of  the  acoustic  wave  speed,  u  +  a,  and  the  waves  convected  at  the  fluid  speed,  u.  The 
application  of  preconditioning  changes  the  eigenvalues  of  the  system  of  compressible  flow  equations 
and  reduces  this  disparity  in  the  wave  speeds.  For  example,  the  time  derivatives  are  premultiplied 
by  a  matrix  that  slows  the  speed  of  the  acoustic  waves  relative  to  the  fluid  speed. 

The  preconditionings  that  are  applied  here  not  only  accelerate  the  convergence  to  a  steady 
state  but  can  also  change  the  steady-state  solution  because  of  the  choice  of  artificial  viscosity,  or 
upwinding,  terms.  Similarly,  the  boundary  conditions  are  based  on  the  preconditioned  equations 
rather  than  the  original  governing  equations.  As  discussed  in  ref.  [9],  the  “standard”  numerical 
schemes  for  the  compressible  equations  do  not  converge  to  the  solution  of  the  incompressible  equa¬ 
tions  (using  a  pseudo-compressibility  approach)  as  the  Mach  number  approaches  zero.  However, 
the  use  of  a  proper  preconditioning  leads  to  a  numerical  scheme  that  does  behave  appropriately  for 
low  Mach  numbers. 

In  this  paper,  we  present  a  generalization  of  the  preconditioners  given  by  Turkel  [10]-  [11],  and 
Choi  and  Merkle  [1],  as  well  as  those  presented  more  recently  by  Radespiel  and  Turkel  [6]  and 
Radespiel  et  al.  [7].  We  discuss  both  nonconservative  and  conservative  artificial  dissipation  models 
and  the  effects  of  the  preconditioning  matrix  on  the  accuracy.  Numerical  results  indicate  that  by 
using  the  conservative  formulation  of  artificial  dissipation  model,  accurate  solutions  are  obtained 
without  sacrificing  efficiency. 

We  show  that  preconditioning  can  be  combined  with  well-known  convergence  acceleration  tech¬ 
niques  such  as  residual  smoothing  and  multigrid.  Indeed,  the  clustering  of  eigenvalues  with  pre¬ 
conditioning  improves  the  damping  of  transient  high-frequency  modes  to  an  extent,  which  makes 
efficient  multigrid  computation  of  low  Mach  number  flows  practical. 


Algorithm 

The  conservation-law  form  of  the  Euler  equations  can  be  transformed  easily  into  non-conservation 
form  by  matrix  transformations  and  vice-versa.  For  convenience,  we  start  with  the  non-conservation 
form  of  the  Euler  equations.  Note  that  although  the  theory  is  developed  for  the  Euler  equations,  the 
methodology  is  applied  in  a  straight-forward  manner  to  the  Navier-Stokes  equations  by  grouping 
the  viscous  fluxes  with  the  dissipative  fluxes. 

We  consider  the  preconditioned  Euler  equations  written  as 

p-'^Qt  +  AQ:c  +  BQy  +  CQz  =  0  (1) 
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The  form  of  the  matrices  and  C  depends  on  the  choice  of  variables  Q.  We  first  consider 

the  variables  Q  =  Qo  =  (p,  u,  v,  w,  S),  where  the  entropy  satisfies  the  relation  dS  =  dp-  a'^dp.  We 
then  have  the  following  form  of  the  matrices  for  Qo- 
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In  generalized  coordinates,  we  are  interested  in  combinations  of  these  matrices.  Hence,  we 
define  Dq  =  Ao^i  +  BqU}2  +  and  q  =  uu:i  +  where  wi,  0:2,  and  ws  are  the  metrics 

associated  with  the  coordinate  transformation.  This  definition  leads  to 
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We  consider  the  preconditioner  Pq  given  by 
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where  a,  13,  and  6  are  free  parameters. 

For  optimal  preconditioning,  should  be  proportional  to  the  square  of  the  local  speed,  u  + 
+  XD^,  ([11]).  However,  this  strategy  introduces  a  complication  near  the  stagnation  points;  the 
preconditioner  becomes  singular  when  /?  =  0.  Furthermore,  Darmofal  and  Schmid  [3]  have  shown 
that  the  eigenvectors  become  less  orthogonal  as  /?  goes  to  zero.  We  introduce  a  simple  cutoff  to  avoid 
this  situation.  Because  this  preconditioner  is  introduced  mainly  for  low-speed  regions,  we  design 
the  preconditioner  to  turn  off  at  higher  speeds.  If  the  preconditioner  is  turned  off  at  a  subsomc 
Mach  number  we  can  use  a  nonconservative  formulation  for  the  artificial  viscosity  without  a  loss  of 
accuracy  in  capturing  weak  solutions.  As  shown  later,  this  strategy  simplifies  the  construction  of 
an  artificial  viscosity.  For  sufficiently  high  Mach  numbers,  we  want  to  remove  the  preconditioning, 
i.e.  0^  =  a^,a  =  0,  and  ^  =  0.  One  choice  is 
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For  nonorthogonal  grids,  +  v'^  +  vP'  can  be  replaced  by  the  sum  of  the  squares  of  the  normalized 
contravariant  velocity  components. 

The  case  of  no  preconditioning  {P  -  I)  corresponds  to  a  =  0,P^  =  a  ,  and  ^  =  0.  The 
parameter  0^  is  returned  to  its  nonpreconditioned  value  at  Mq,  which  is  the  cutoff  value  for  the 
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Mach  number.  Numerical  evidence  suggests  that  K2  depends  on  the  number  of  mesh  nodes  near 
the  stagnation  point  and  possibly  should  depend  on  the  local  cell  Reynolds  number  [1].  Typically, 
Ki  is  between  1  and  1.1,  and  is  between  0.4  and  1.  The  eigenvalues  of  the  matrix  PqDq  are 

Ao  =  q,q,q  (repeated  eigenvalues) 


and 
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So  A  =  dia5(A+,A_,Ao,  Ao,  Ao). 

The  eigenvalues  of  Fo-Do  are  independent  of  6;  however,  the  same  is  not  true  for  the  eigenvectors. 
Though  Choi  and  Merkle  [1]  employ  ^  =  1,  we  consider  ^  =  0  because  this  simplifies  the  eigenvectors 
of  PoDq.  The  right  eigenvectors  are  given  by  the  columns  of 
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Note  that  (A+  -  q){X-  -  q)  =  -  (3'^\u)\‘^.  We  define  some  temporary  quantities.  Let 
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Note  that  R  is  singular  for  =  0,  which  is  merely  an  artifact  because  multiple  eigenvalues  exist 
in  a  multidimensional  eigenspace.  If  U2  =  0,  then  we  can  change  the  eigenvectors  in  the  invariant 
subspace  of  the  multiple  eigenvalues  so  that  either  wi  or  013  appears  in  the  denominator.  Because 
aU  three  of  these  eigenvalues  cannot  be  zero  simultaneously,  some  set  of  nonsingular  eigenvectors 
always  exists.  Also  note  that  UJ2  does  not  appear  in  the  denominator  in  the  final  analysis  and  does 
not  create  any  numerical  difficulties;  therefore  we  can  ignore  this  anomaly. 

The  largest  eigenvalue  of  PqDq  is  used  to  determine  the  (inviscid)  time  step,  ^or  M  ~  0, 
A±  ((1  -  a)  ±  v'(l  -  0)2  +  4^  For  a  =  0,  the  condition  number  of  PqDq  is  ~  2.6,  and 

for  a  =  1  the  condition  number  is  1. 

The  above  matrices  were  given  for  Qo  =  (p, «,  v,  w,  S)  variables.  In  the  code  we  base  everything 
on  Qa  =  (p,  u,  V,  w,  T)  variables.  Then, 

p  _dQApdQo 

dQo  °dQ4 

The  transformations  between  the  Qq  and  Qa  variables  are  given  in  the  appendix. 


Artificial  Viscosity 


For  a  central- difference  scheme,  it  is  necessary  to  add  artificial  dissipative  terms  to  the  finite- 
difference  approximation  of  the  spatial  derivatives  to  damp  the  numerical  oscillations.  A  nonlinear 
second-difference  term  is  normally  added  to  control  oscillations  near  shocks,  and  a  linear  fourth- 
difference  term  is  added  to  damp  high-frequency  oscillations  [5].  We  are  interested  primarily  in 
the  functional  form  of  these  differences.  Hence,  our  examples  include  a  second-difference  artificial 
dissipation;  extensions  to  fourth  differences  and  nonlinearities  are  straight-forward.  Similarly,  for 
one-sided  schemes  the  central  difference  plus  the  artificial  dissipation  is  replaced  by  a  Roe  matrix 
formulation. 

A  typical  preconditioned  finite-difference  scheme  can  be  expressed  as 


AQc  =  AtPc 
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where  Qc  represents  the  conservative  variables  and  the  spatial  partial  derivatives  are  replaced  by 
central-difference  approximations.  We  consider  both  conservative  and  non- conservative  ways  of 
adding  artificial  dissipative  terms.  The  dissipation  need  not  be  expressed  in  terms  of  the  conserva¬ 
tive  variables  Qc-  For  another  set  of  variables  Qv,  we  get  AQc  =  §^AQv-  We  choose  our  “basic” 
form  as  the  one  given  in  Qo  =  {p,u,v,w,S)  variables  (i-e-,  Fq).  Then,  for  another  set  of  variables 

Qv,  we  have  the  preconditioner  Fy  =  =  a^Fog^. 

We  add  second-difference  terms  for  the  artificial  viscosity  to  Eq.  (6)  and  omit  the  fourth 

differences  for  brevity.  We  now  get 
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where  Fy  =  Fy^^.  Thus,  Fy  contains  both  the  preconditioner  and  the  change  of  variables.  The 
above  formulation  is  non- conservative  because  Fy  is  outside  the  derivative  terms.  In  the  above 
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equation,  c  is  a  matrix  function.  For  example,  if  a  is  proportional. to  the  spectral  radius,  then  we 
have  a  scalar  artificial  viscosity  in  Qv  variables;  whereas  if  cr{A)  ~  A,  then  we  have  a  matrix- valued 
artificial  viscosity. 

After  we  have  computed  AQy,  we  perform  the  residual  smoothing  on  it.  In  the  multigrid 
algorithm,  AQv  is  passed  to  the  next  coarser  grid.  At  the  end  of  each  stage,  Qv  is  recomputed. 
The  conservative  variables  are  then  calculated  as  nonlinear  functions  of  the  Qv  variables.  We 
should  mention  that  Av  —  Ac  i  similar  expressions  can  be  written  for  the  other  coordinate 

directions.  Let  Av  —  qq^  be  the  Jacobian  matrix;  then,  Av  =  qq^-^v-  Hence, 


FvAv  = 


dQv p  .  dQc 
dQc  ^  ‘^dQv 


=  FyAc 


dQc 

dQv 


Because  the  spectral  radius  is  invariant  under  a  similarity  transformation  for  a  scalar  viscosity, 
we  can  use  either  PvAv  or  PcAc-  Let  A  =  diag(Ax,  A2,A3,A3,  A3),  where  Ai,A2,  and  A3  represent 
A+,  A_,  and  Aq  (i.e.,  the  eigenvalues  of  PD),  modified  by  cutoffs  near  the  stagnation  points.  Define 
|A|  =  diag(|Ai|,|A2|,|A3|,|A3|,|A3|).  Then,  |PD|  =  R-^\A\R.  For  any  vector  a;  =  {xi,X2,X3,X4,X5y, 
\PD\x  =  (P"^|A|)z,  where  z  =  Rx.  Define  Y  =  u;iX2  -f  0J2X3  +  ^3X4.  Then, 
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where  zi  is  the  first  element  of  the  vector  z.  We  note  that  Z3  and  Z4  do  not  appear  by  themselves 
in  subsequent  formulas  but  rather  in  conjunction  with  other  variables  in  the  form 
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For  a  scalar  viscosity  (i.e.,  when  a  is  the  spectral  radius),  the  viscosity  remains  scalar  after 
preconditioning.  In  this  approach,  differences  of  the  Qv  variables,  rather  than  those  of  the  conserved 
variables,  are  added  to  each  of  the  conservation  equations. 
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However,  the  formulation  used  until  now  is  nonconservative.  Since  we  are  also  interested  in 
solving  transonic  flows  with  shocks  using  this  scheme,  we  present  a  conservative  formulation  of  the 
artificial  dissipation.  By  dropping  the  subscript  V ,  we  obtain 


AQ^AtT 


d£  ^  dH 
dx'^  dy  ^  dz 


X 


Even  if  (7  is  a  scalar  function,  we  must  evaluate  a  matrix- vector  product;  as  a  result,  the  numerical 
effort  is  equivalent  in  complexity  to  a  matrix- valued  viscosity. 

Following  previous  work  on  multigrid  schemes  for  the  Navier- Stokes  equations,  these  artificial 
viscosity  functions  account  for  the  ratio  of  the  spectral  radii  in  the  different  coordinate  directions, 
that  is. 


a{PA)  =  k^''MP^) 


(p{PB)\^  (p{PC)\^' 

\p{PA)j  ^\p{PA)). 


where  d  denotes  the  original  viscosity  function,  p  denotes  the  spectral  radius,  and  ^  is  the 
artificial  viscosity  coefficient  that  corresponds  to  the  fourth-difference  terms.  Similar  relations  are 
used  to  modify  the  artificial  viscosity  by  cell  aspect-ratio  in  the  other  coordinate  directions.  These 
scaling  functions  provide  sufficient  artificial  dissipation  for  general  curvilinear  grids  that  contain 
high  aspect-ratio  cells. 

In  the  previous  section,  a  preconditioner  was  introduced  that  is  dependent  on  the  parameters 
a  and  6.  Because  6  does  not  affect  the  eigenvalues  of  PA,  it  has  no  affect  on  the  scalar  artificial 
viscosity.  Choosing  a  =  1  reduces  the  largest  eigenvalue,  which  also  improves  the  condition  number 
and  decreases  the  artificial  viscosity  compared  with  a  =  0  case.  We  thus  expect  a  =  1  to  slow  the 
convergence  compared  with  the  case  in  which  a  =  0,  in  spite  of  the  fact  that  we  have  reduced  the 
condition  number.  However,  we  expect  that  with  a  =  1,  the  numerical  accuracy  will  improve.  For 
a  matrix  viscosity  (or  a  Roe  matrix),  we  expect  similar  but  less  pronounced  dependence  on  o:. 


Boundary  Conditions 

In  many  CFD  codes,  the  boundary  conditions  in  the  far  field  are  based  on  characteristic  variables, 
even  for  viscous  flow.  Thus,  at  inflow  the  incoming  variables  that  correspond  to  positive  eigenvalues 
are  specified,  and  the  outgoing  variables  that  correspond  to  negative  eigenvalues  are  extrapolated. 
A  change  in  the  time-dependent  equations  also  changes  .the  characteristics  of  the  system  (although 
the  signs  of  the  eigenvalues  remain  unchanged).  Hence,  the  boundary  conditions  must  be  modified 
for  the  preconditioned  system. 

In  the  present  study,  we  have  used  the  simplified  far-field  boundary  conditions  suggested  in  ref. 
[6].  Basically,  all  variables  at  far-field  boundaries  are  specified  in  terms  of  two  sets  for  p,u,v,w, 
and  T.  Depending  on  whether  the  subsonic  flow  is  an  inflow  or  an  outflow,  one  set  is  specified 
at  the  free-stream  levels,  and  the  other  set  is  extrapolated  from  the  interior.  For  example,  at  the 
inflow 

m  =  Woo?  ^6  =  Wb  =  Woo?  Tt  =  TooyVh  =  Pint 
and  at  the  outflow 

Uh  =  Uint^  Vh  =  Vint^'^h  =  ^intj  =  Tini^Ph  =  Poo 
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where  the  subscripts  “6”  and  “inf’  refer  to  the  values  at  the  boundary  and  adjacent  interior  points, 
respectively. 

For  supersonic  flows,  standard  boundary  conditions  are  used;  that  is,  extrapolation  at  the 
outflow  boundary  and  free-stream  values  for  all  variables  at  the  inflow  boundary  are  prescribed. 


Changes  to  Original  Coding 

Here,  we  present  the  steps  necessary  to  introduce  preconditioning  into  an  existing  compressible  flow 
code.  We  assume  an  explicit  time-stepping  scheme  (e.g.,  a  Runge-Kutta  scheme)  that  is  augmented 
by  an  implicit  residual  smoothing  and  multigrid  scheme  to  obtain  steady-state  solutions. 

1.  Precondition  the  residual. 

•  Select  a  reasonable  set  of  dependent  variables  for  the  problem  (e.g.,  p,  u,u,w,T).  Note 
that  other  sets  of  variables  are  possible. 

•  For  nonconservative  scalar  artificial  viscosity 

-  Multiply  the  physical  residual  (inviscid  and  Navier-Stokes  portions  and  physical 
forcing  functions  (if  applicable))  by  F  =  Pt^^- 

-  Add  scalar  artificial  viscosity  in  these  new  variables  to  the  residual. 

•  For  conservative  scalar  viscosity 

-  Add  inviscid,  viscous,  and  artificial  viscosity  fluxes  to  obtain  the  total  residual. 
Scalar  viscosity  already  includes  the  transformation  F”^ . 

-  Multiply  this  total  residual  by  F. 

•  Multiply  the  total  preconditioned  residual  in  the  (p,  u,  v,  w,  T)  by  At.  Apply  the  residual 
smoothing  for  the  changes  in  terms  of  the  (p,u,v,w,T)  variables.  Add  the  residuals  in 
the  (p,  u,  V,  w,  T)  variables  to  the  dependent  variables  at  the  previous  time  iteration 
to  obtain  the  values  of  {p,u,v,w,T)  at  the  next  time  step.  The  new  values  of  the 
conservative  variables  are  then  evaluated  as  nonlinear  functions  of  (p,  u,v,w,T). 

2.  For  the  FAS  multigrid  algorithm,  we  need  to  restrict  the  residuals  and  the  variables  from  a 
finer  grid  to  a  coarser  grid.  In  the  computer  code,  the  residuals  are  stored  in  the  (p,  u,  u,  w,  T) 
frame,  whereas  the  dependent  variables  are  stored  as  conservative  variables.  Hence,  restric¬ 
tions  and  prolongations  are  done  on  the  conservative  variables.  The  residual  smoothing  is 
performed  on  the  (p, «,  v,  w,  T)  residuals  and  the  forcing  functions  on  the  coarse  grids  are 
evaluated  in  (p, «,  v,  w,  T)  variables. 

.3.  Choose  a  new  time  step  for  the  inviscid  portion  based  on  the  preconditioned  system.  The 
viscous  contribution  to  the  time  step  is  then  incorporated  as  usual. 

4.  Modify  the  far-field  boundary  conditions. 

Computational  Results 

An  existing  three-dimensional  compressible  Navier-Stokes  flow  solver  was  modified  to  include  the 
preconditioning  methodology  described  in  the  preceding  paragraphs.  The  modified  code  was  used 
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Mach  no. 

panel  method 

no  preconditioning 

preconditioning 

Cl 

0.1 

.241 

.2448 

.24.34 

0.01 

.2172 

.2421 

0.001 

.1097 

.2421 

Cd 

0.1 

0.0 

.0008 

.00027 

0.01 

.0086 

.00027 

0.001 

.0731 

.00027 

Table  1:  Comparison  of  lift  and  drag  with/without  preconditioning 

to  compute  low-speed  and  transonic  flow  over  several  configurations  of  practical  interest.  Three  of 
these  test  cases  are  discussed  here.  The  computational  results  reported  in  this  paper  were  obtained 
with  scalar  form  of  artificial  dissipation. 

Two  Dimensional  Results 

We  first  consider  a  two  dimensional  inviscid  subsonic  problem.  For  this  case  we  can  compare  the 
compressible  code  with  and  without  preconditioning  to  a  panel  method  in  order  to  assess  their 
accuracy.  We  consider  inviscid  flow  over  NACA  0012  airfoil.  The  C-type  grid  has  224  x  40  cells 
clustered  near  the  leading  and  trailing  edges  in  order  to  allow  accurate  drag  computations.  As  seen 
in  table  1  the  code  using  the  preconditioning  gives  lift  and  drag  quite  close  to  the  panel  method 
results  with  only  a  small  dependence  on  the  Mach  number.  The  non-preconditioned  code  has  large 
variations  as  the  Mach  number  goes  to  zero  and  is  not  converging  to  the  panel  code  results.  In 
particular  the  drag,  which  should  be  zero  for  inviscid  flow,  is  very  large  without  preconditiomng 
but  close  to  zero  when  using  preconditioning. 

Low-Speed  Flow  over  Three-Dimensional  Wing 

Essentially  incompressible  viscous  flow  over  the  ONERA  M6  wing  is  considered  as  the  first  test 
case.  For  this  case,  the  Reynolds  number  (based  on  mean  aerodynamic  chord)  is  11.7  million, 
and  the  angle  of  attack  is  .3.06°.  A  grid  that  consists  of  193  x  49  x  33  points  is  used  for  these 
computations. 

The  first  set  of  results  in  Fig.  1  shows  the  effect  of  preconditioning  on  the  convergence  history 
of  the  numericcd  algorithm  at  a  Mach  number  of  0.1.  For  these  computations,  50  iterations  on  the 
coarse  grid  are  Mowed  by  300  iterations  on  the  fine  grid.  This  figure  clearly  shows  a  significant 
improvement  in  the  convergence  rate  when  preconditioning  is  used. 

The  effect  of  the  free-stream  Mach  number  is  considered  in  the  next  series.  Figures  2  -  3 
show  the  effect  of  free-stream  Mach  number  on  convergence  rates  and  surface  pressure  distribution, 
respectively.  The  residuals  in  Fig.  2  have  been  normalized  with  their  respective  initial  values  to 
remove  the  scaling  effects  caused  by  differences  in  the  free-stream  Mach  number.  Note  that  over 
a  Mach  number  range  of  0.01  to  0.2,  the  convergence  rates  for  the  preconditioned  scheme  are  very 
similar;  the  asymptotic  convergence  rates  for  the  two  lowest  Mach  numbers  are  almost  identical. 
Although  not  shown  here,  similar  results  have  been  obtained  at  an  even  lower  Mach  number  of 


8 


0.001.  On  the  other  hand,  the  original  non-preconditioned  scheme  failed  to  converge  at  Mach 
numbers  of  0.01  and  lower;  similar  to  the  observations  reported  by  Volpe  [12]. 

The  pressure  distributions  (at  a  span  location  of  80%)  shown  in  Fig.  3  indicate  that  results  for 
Mach  numbers  of  0.01  and  0.1  are  identical  within  plotting  accuracy.  Even  at  a  Mach  number  of 
0.2  (except  for  a  slightly  higher  value  of  the  pressure  peak  in  the  leading-edge  region),  the  effect  of 
Mach  number  is  negligible.  These  results  demonstrate  that  the  preconditioned  system  approaches 
the  incompressible  limit  in  a  smooth  and  systematic  manner  without  any  penalty  in  convergence 
rate. 

The  effects  of  the  conservative  and  non- conservative  form  of  the  artificial  viscosity  formulations 
are  shown  in  Figs.  4  and  5.  The  non-conservative  formulation  converges  better  on  the  coarser  grid, 
but  is  nearly  identical  in  performance  to  the  conservative  formulation  on  the  fine  grid.  The  effect 
of  these  two  forms  of  artificial  viscosity  on  pressure  distribution  is  negligible,  which  is  expected  to 
be  the  case  at  low  speeds. 

Flow  over  Two-Dimensional  Multi-element  Airfoil 

The  next  test  case  considered  here  is  that  of  a  3-element  airfoil  configuration  that  has  been  inves¬ 
tigated  both  experimentally  and  theoretically  [13].  Present  computations  are  performed  at  a  chord 
Reynolds  number  of  9  million  and  an  angle  of  attack  of  16.2°.  The  Mach  number  for  this  test  case 
is  0.2.  A  20-block  structured  grid,  shown  in  Fig.  6,  is  used  for  the  computations.  This  test  case 
has  also  been  investigated  with  the  original  non-preconditioned  version  of  the  flow  solver  used  in 
this  work,  as  reported  by  Vatsa  et  al.  [14] . 

The  convergence  histories  for  this  case  are  presented  in  Fig.  7,  where  the  results  from  the 
original  non-preconditioned  and  current  preconditioned  (conservative)  schemes  are  compared.  The 
residuals  in  the  original  scheme  indicate  considerable  slowdown  in  convergence  at  approximately  4- 
orders;  whereas  the  residuals  in  the  preconditioned  scheme  exhibit  much  better  convergence.  Note 
that  the  free-stream  Mach  number  of  0.2  for  this  case  is  not  considered  too  low  for  compressible 
codes.  However,  in  this  case  several  pockets  of  slow-moving  flow  exist  in  the  cove  regions  of  the 
slat  and  the  main  airfoil  sections  [14];  these  pockets  slow  the  convergence  of  standard  compressible 
codes.  The  preconditioned  system  has  a  much  better  condition  number  for  the  eigenvalues  and 
does  not  experience  slowdown  as  a  result  of  such  disparities  in  the  flow  speed. 

The  computed  pressure  distributions  for  this  case  are  compared  in  Fig.  8.  As  expected,  little 
difference  is  observed  in  the  two  sets  of  computed  results;  furthermore,  these  results  compare  quite 
favorably  with  the  experimental  data. 

Transonic  Flow  over  Three-Dimensional  Wing 

The  final  test  case  presented  here  involves  transonic  flow  over  the  ONERA  M6  wing.  The  test 
conditions  for  this  case  are  identical  to  the  first  test  case  except  for  the  free-stream  Mach  number, 
which  is  chosen  as  0.84.  The  Navier-Stokes  solutions  for  this  case  were  obtained  for  conservative  and 
non- conservative  preconditioners.  A  baseline  solution  with  no  preconditioning  was  also  obtained 
for  comparison. 

The  computed  pressure  distributions  for  this  case  are  compared  in  Fig.  9  at  80%  span  station. 
This  figure  clearly  shows  that  the  pressure  distributions  from  conservative  preconditioning  are  vir¬ 
tually  indistinguishable  from  the  baseline  (unpreconditioned)  case.  However,  the  non-conservative 
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form  of  preconditioning  produces  noticeable  differences  in  the  pressure  distribution  in  the  vicinity 

of  the  shocks.  .  .  ,  j  v  *  *1.  v 

The  convergence  histories  for  this  case  are  shown  in  Fig.  10,  where  it  is  observed  that  the  base¬ 
line  scheme  stalls  once  the  residual  drops  approximately  4-orders  in  magnitude.  The  convergence 
properties  of  the  two  preconditioned  schemes  are  similar  to  one  another  and  better  than  those  of 
the  baseline  scheme.  Note  that  the  slow  convergence  of  the  basic  scheme  cannot  be  attributed  to 
low  free-stream  velocities.  However,  the  computational  domain  consists  of  low-speed  flows  in  the 
stagnation  and  boundary-layer  regions,  where  the  convective  speed  of  propagation  is  much  slower. 
Preconditioning  appears  to  reduce  the  imbalance  of  convective  and  acoustic  speeds  in  these  regions, 
which  improves  the  overall  convergence  rate  for  such  problems. 

Concluding  Remarks 

An  attractive  scheme  for  computing  low-speed  flows  has  been  presented  here  in  the  framework  of 
a  preconditioning  applied  to  the  compressible  flow  equations.  The  current  formulation  produces 
accurate  incompressible  results  in  a  smooth  and  systematic  manner  as  the  Mach  number  approaches 
zero.  Moreover,  the  proposed  scheme  is  relatively  easy  to  implement  in  an  existing  compressible 

flow  code.  j-x-  j 

The  test  cases  presented  here  demonstrate  the  efficiency  and  accuracy  of  the  preconditioned 

scheme.  Excellent  convergence  has  been  obtained  at  Mach  numbers  that  range  from  0.01  to  0.84. 
The  resulting  pressure  distributions  agree  well  with  known  solutions  for  these  cases.  Based  on  our 
experience  thus  far,  this  scheme  offers  a  viable  alternative  to  purely  incompressible  flow  codes  for 
computing  low-speed  flows.  In  addition,  this  scheme  offers  the  advantage  of  being  able  to  compute 
flows  with  mixed  speed  regimes,  in  which  the  local  Mach  numbers  Can  vary  from  very  low  subsomc 
to  supersonic  values,  e.g.,  in  the  case  of  flow  over  high-Hft  configuration  near  maximum  lift.  FinaUy, 
this  scheme  can  improve  the  convergence  rate  even  for  viscous  transonic  flows  by  preconditioning  the 
embedded  low-speed  flows  in  the  boundary  layers.  Future  work  should  focus  on  matrix-dissipation 
or  Roe  type  schemes,  to  improve  the  numerical  accuracy. 


Appendix 

We  define  the  sets  of  variables 


Qo  =  {p,u,v,w,S) 


Qc  =  Qi  =  {p,pu,pv,pw,E) 


and 


Q4  =  (p,u,v,w,T). 

Let  a  be  the  speed  of  sound,  and  q'^  =  u'^  +  v^  +  w'^.  Given  the  preconditioning  Pq  in  Qo 
variables,  we  can  compute  the  preconditioner  Pi  in  Qi  coordinates  by 


,  dQi  dQo 
dQo  ""dQi^ 


where 


Po  is  given  by  Eq.  (2). 
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The  following  transformation  matrices  connect  these  variables;  Let  ^  y^-- 


dQo 

dQi 


—u 

p 

p 

’-vu 


^  0 


(l-7)w  (1-7)1;  (1— 7)u;  7-I 

i  0  0  0 

0  T  0  0 

(l-7)w  (1—7)1;  (1— 7)n;  7-I  J 


1 

p 

0 


dQo 


1 

u 

a? 

V 

w 

n2 


0 

p 

0 

0 


0 

0 

p 

0 


1  I 

T-1  2 


0 

0 

0 

p 


pu  pv  pw 


—u 

—  V 


2  J 


dQo 

QQa 


1  0  0  0 
0  10  0 

0  0  10 

0  0  0,1 

1-7  000 


0 

0 

0 

0 

TP 

T 


dQA 

dQo 


1 
0 
0 

0  0  0  1 

0  0  0 

TP 


0 

0 


0  0  0 
1  0  0 
0  10  0 
0 


n 

TP 


To  transform  the  preconditioned  residual  in  (p,  u,  u,  w,  T)  to  and  from  the  conservative  variables, 
we  use  the  following  Jacobians,  which  result  directly  from  those  given  above: 


dQ\ 


p 

p 

0 

0  0 

-P 

T 

dQ^_ 

pu 

p 

pv 

p 

0 

0 

0  0 

P  0 

0  P 

—pu 

T 

—pv 

BQa 

p 

pw 

p 

T 

—pw 

T 

E 

_  P 

pu 

pv  pw 

-pq^ 

2T 

(1- 

7)w 

(l-7)t; 

(1-7) 

_ u 

i 

0 

0 

0 

0 


1 

P 

0  • 
p 


0 

i 

p 

P 


0 

0 

0 

p 


Also,  Ti  =  so  that  Pi  =  ri|§f.  In  particular,  Pc  =  |§7r4  and  P4  =  r4^.  For  the 


'dQi 


dQP 
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{p,u,v,  w,T)  variables,  we  consider  only  the  case  in  which  ^  =  1. 


We  define  f  =  Then, 


r4 


0^ 

^(l  +  «) 

^(l  +  o) 


0  0  0  0' 

1  0  0  0 

0  i  0  0 

0  0  i  0 

-fu  -fv  -fw  f ' 


and  letting  E  =  : 


0  0  0  0  ■ 

p  0  0  0 

0  /)  0  0 

0  0  P  0 

j^E  pu  pv  pw  _ 
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Figure  1:  Effect  of  preconditioning  on  convergence  history  for  ONERA  M6  wing  at  Moo  -  0.10 


Figure  2:  Effect  of  Mach  number  on  convergence  history  of  the  preconditioned  scheme  for  ONERA 
M6  wing 
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-1.5 


- M=0.20,  cons.  prec. 


Figure  3:  Effect  of  Mach  number  on  ONERA  M6  wing  pressure  distributions  at  80%  span  location 
for  preconditioned  scheme 


Figure  4:  Effect  of  conservative/non-conservative  formulations  of  artificial  dissipation  on  conver¬ 
gence  history  for  ONERA  M6  wing 
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- M=0.01,cons.  prec. 

- M=0.01,  non-cons.  prec. 
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Figure  5:  Effect  of  conservative/non-conservative  formulations  of  artificial  dissipation  on  ONERA 
M6  wing  pressure  distributions 


Figure  6:  Partial  view  of  block-structured  grid  for  3-element  high  lift  airfoil  configuration 
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Figure  7:  Effect  of  preconditioning  on  convergence  history  for  the  the  3-element  high-lift  airfoil 
configuration 


Figure  8:  Effect  of  preconditioning  on  pressure  distributions  for  the  3-element  high-lift  airfoil 
configuration 
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x/c 


Figure  9:  Comparison  of  pressure  distribution  for  ONERA  M6  wing  at  transonic  speeds 


Figure  10:  Comparison  of  convergence  histories  for  ONERA  M6  wing  at  transonic  speeds 
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